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Abstract This paper is devoted to the treatment of specific numerical problems which 
appear when phase appearance or disappearance occurs in models of two-phase flows. 
Such models have crucial importance in many industrial areas such as nuclear power 
plant safety studies. In this paper, two outstanding problems are identified: first, the 
loss of hyperbolicity of the system when a phase appears or disappears and second, the 
lack of positivity of standard shock capturing schemes such as the Roe scheme. After 
an asymptotic study of the model, this paper proposes accurate and robust numerical 
methods adapted to the simulation of phase appearance or disappearance. Polynomial 
solvers are developed to avoid the use of eigenvectors which are needed in usual shock 
capturing schemes, and a method based on an adaptive numerical diffusion is designed 
to treat the positivity problems. An alternate method, based on the use of the hyper- 
bolic tangent function instead of a polynomial, is also considered. Numerical results are 
presented which demonstrate the efficiency of the proposed solutions. 

Key words: two-phase flows, numerical simulation. Roe scheme, hyperbolic system, 
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1 Introduction 



Multiphase flows can be found in a large variety of industrial or natural systems involving 
boiling or condensing fluids, reacting flows or aerosols. Such systems are, e.g., power 
plants, refrigerators, distillation units, gas or oil pipelines, pollutant separators, or clouds. 
The present work has been conducted in the context of nuclear power plant safety studies. 
In nuclear reactors, the appearance of vapor around the fuel rods interferes with the heat 
evacuation and can cause severe damages. To design and optimize the equipments in 
order to guarantee the highest possible safety level, numerical simulations of multiphase 
flows are intensively used. However, these simulations remain extremely delicate because 
of the complexity of the models and the possible huge discrepancy between the volume 
fraction of the various phases. For instance, within a subcooled liquid injected in a heated 
column, a transition from a single-phase liquid at the inlet to pure vapor at the outlet 
may take place. In such situations, numerical difficulties may be observed, like the loss of 
positivity of the mass fractions or internal energies. This is the case for instance with the 
CEA research code OVAP [21] based on an implicit version of Roe's scheme. Therefore, 
a robust numerical scheme for two-phase flows must be able to treat all ranges of volume 
fractions. 

In the literature, few works deal with the problem of phase appearance and disap- 
pearance explicitly. In most codes, this problem is treated using ad-hoc fixes. A first 
treatment has been developed in the "CATHARE" two-phase fiow code [3]. It relies on 
specific expressions of the interfacial mass and energy transfer terms which are designed 
such that the void fraction remains in an interval [oimin-, otmax]- This treatment is com- 
bined with a numerical conditioning of the interfacial and wall friction source terms in 
order to provide a proper mechanical model for the coupling of the residual phases. A 
similar strategy is used in the "NEPTUNE" code [17]. A second method is proposed in 
|29j where an extension of the AUSM+ scheme (Advection Upwind Splitting Methods) 
to two-fiuid models is developed. In [29], it is noticed that the AUSM+ numerical fiuxes 
remain non-singular when a phase disappears as long as the involved Mach number and 
phase velocities remain bounded. Since it is assumed that the velocity of the two phases 
should tend to each other at the transition, the velocities are therefore artificially tied to 
each other through a smooth function. A similar treatment is applied to the temperature. 
This treatment is applied when a phase has a volume fraction below amin = 10~^. 

Therefore, the strategy developed in the literature is to treat these problems at the 
level of the underlying physics, by designing specific expressions for the interfacial closure 
terms. Without underestimating the role of the physics, we propose an alternative route. 
We explore the mathematical structure of the two-phase models in the limit of small 
volume fraction of one the phases. This asymptotic approach is used to highlight the 
possible causes of the numerical breakdown and to design more robust methods. We 
restrict ourselves to models of two-phase fiows but the methods could be extended to 
three of more phases. In the numerical investigations, we rely on a time-implicit version 
of the Roe scheme used in the "OVAP" code [24J. We identify two essential difficulties, 
which are: (i) the loss of hyperbolicity of the two-fiuid model when a phase appears or 
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disappears, and (ii): the lack of positivity of the Roe scheme. Each of these difficulties 
will receive a specific treatment. 

To address the first difficulty, we propose the use of the so-called polynomial schemes 
|10j . This choice is motivated by the asymptotic analysis of the two-phase model in the 
limit of vanishing volume fraction of one of the phases. In this limit, the two phases almost 
decouple and the minority phase obeys a pressureless gas dynamics system [HI E] . This 
system is not hyperbolic because the Jacobian of the fiux matrix is not diagonalizable. 
This implies that two eigenvectors of the original two-phase model collapse in the limit of 
small volume fraction. Therefore, most shock-capturing schemes, which require a complete 
basis of eigenvectors, breakdown in this limit. To overcome this problem, schemes that 
do not require that eigenvectors form a complete basis are needed. There are many such 
schemes, such as Lax-Friedrichs, or central schemes [2B], but many of them are too diffusive 
for safety studies of nuclear power plants. The interesting feature with the polynomial 
schemes is that it is possible to tune the amount of numerical diffusion. Polynomial 
schemes have been used e.g. in [2U 128] . 

We will also consider an alternate method which uses the hyperbolic tangent function 
instead of a polynomial. It is as precise as the most precise polynomial method, and shows 
very good positivity properties without requiring any positivity treatment (see below). 
However, it is currently computationally too intensive for practical use. Nonetheless, 
improvements in the efficiency of the computation of the hyperbolic tangent function of 
a matrix could make this method potentially very competitive. 

The second difficulty, namely, the lack of positivity, is a critical issue in phase-transition 
problems. Indeed, they frequently appear in areas where the mass fraction of one of 
the phases is small. Then, small inaccuracies easily lead to negative mass fractions, 
especially with large time-steps. The simple fix consisting in replacing negative quantities 
by arbitrarily small positive values results in conservation losses and degraded robustness. 
For this reason, the development of positive schemes has been considered a major issue. 
In [2], Einfeldt and al introduced the notion of "positively conservative" schemes where 
the density and internal energy remain positive. While the Godunov scheme is positively 
conservative, they show that no linearized Riemann solver, including the Roe scheme, 
is positively conservative. A more detailed bibliography about positive schemes can be 
found in section IH 

Several specific aspects make previously developed methods of difficult use for two- 
phase fiow models, especially in the context of nuclear power plants safety. First, the 
models, such as the ones presented in the forthcoming sections, are complex non con- 
servative systems. The analytical expressions of the eigenvalues and eigenvectors are not 
available. An analytical proof of the positivity of a scheme is therefore not possible. Addi- 
tionally, such proofs strongly use the eigenstructure of the system. However, as explained 
above, this eigenstructure becomes singular when the volume fraction of one of the phases 
becomes small. As we will see in the review of section |4[ many strategies leading to pos- 
itive schemes are based on an increase of the numerical diffusion. But this additional 
diffusion is detrimental for the accuracy of the scheme, and accuracy is a critical issue for 
the targeted application. Another critical issue is efficiency and motivates the use of im- 
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plicit schemes and large time-steps. In this context, schemes inducing positivity through 
a restriction on the CFL stabihty condition are not acceptable either. 

Our method does not guarantee positivity in all cases but, in practice, it solves most 
of the positivity problems while meeting the constraints listed above. The numerical 
treatment consists in an adaptive diffusion, which corrects positivity problems where they 
occur, locally in space and time. It is inspired from the works of Gallice [15] and Romate 
|32j . but the proposed strategy, which uses the framework of the polynomial schemes, 
with a specific choice of the polynomial, is, to our knowledge, original. 

The paper is organized as follows: section [2] develops the asymptotic study of the 
two-phase model when one of the phases disappears, showing that the model loses its 
hyperbolicity in this limit. Section |3] proposes the use of polynomial schemes in replace- 
ment of the Roe scheme to overcome the problem highlighted in section [2j It provides a 
comparison between various choices for the polynomial and selects the most robust one. 
A method similar to polynomial solvers and using the hyperbolic tangent function is also 
detailed. Section |4] addresses the positivity problem and proposes a new strategy to deal 
with it. Finally, numerical results are presented in section [5j A conclusion is given in 
section |6j Auxiliary calculations are collected in appendix \K\ 



2 Two-phase flow models and phase appearance or 
disappearance 

2.1 The full two-phase model 

This paper is concerned with two-phase flow models. Detailed derivations and descriptions 
of two-phase flow models can be found in [191 EO]- In this section, we present the full 
two-phase model, including energy equations, which will be used in the numerical tests of 
section |5| Below, in section 2^, an asymptotic analysis of the simpler, isentropic version 



of this system will be conducted. 

The unknown physical quantities are the volume fraction G [0,1], the density 
Pk > 0, the velocity Uk G M.'^, the energy Ek > 0, the enthalpy hk > of each of the 
phases indexed by k, where the subscript k stands for i for the liquid and v for the vapor. 
They depend on position a; G M*^ (where d is the dimension), and time t. The common 
pressure of the two phases is denoted by p. Here, pressure equilibrium between the two 
phases is postulated. This hypothesis is known as the hydrostatic assumption. The model 
is written as follows, ignoring the viscous terms for simplicity: 

dt{a^Pv) + V ■ {ayp^u^) = T , (2.1) 

dtiaipi) + V ■ (aiPiUi) = -T , (2.2) 
dt{a^PvUy) + V ■ {a^PvUv ® Uy) + a^Vp + DpiVa^ = 

= Tu' + aypJext + i^f + , (2.3) 
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dt{aiPiUi) + V ■ {aepeue ® ue) + agVp + DpiVae = 

= -Tv' + aepife.t + + Fl , (2.4) 
P 

dt{a^PvEy) +pdta^ + V ■ {a^pyUy{Ey H )) = 

Pv 

= T{^ul + hi) + ayPyUt ■Uy + Q: + ■ , (2.5) 

p 

dt{aePiEi) + pdtae + V ■ {aipeUe{Ee H )) = 

Pe 

= -r(^n2 + hi) + aepde.t ■ u, + Q"^, + F^^ ■ u' , (2.6) 
Ot, + = 1 , (2.7) 

2 

Pv = Pv{P, hy), hy=Ey-^-\ , (2.8) 

Z Py 

2 

Pi = pi{p, hi), he = Ei - ^ + ^. (2.9) 

^ Pe 

where Dpi is the interfacial pressure default proposed by Bestion ^ and given by: 

D'p = 6ayaep\ur\^. (2.10) 

The average density p and the relative velocity Ur are defined by 

PvPe 

p = ■ , Ur = Uy- Ue, 

aypi + aipy 

and 5 is an ad-hoc coefficient. Pvip, hy) and pi{p, hi) are the vapor and liquid equations- 



of-state. In the isentropic case (i.e. py = py{p) and pi = pi{p) only, see section 2.2) 



expression (2.10) guarantees that the system is hyperbolic provided that 5 > 1 [36] . 

The source terms have complex physical interpretations and we refer to [20] for 
details. They will not be discussed here. Specifically, F is the interfacial mass transfer 
term, is the interfacial velocity, fext is an external force such as gravity, FJ.^ is the drag 
force, F^ is the wall friction for each phase, h], is the interfacial liquid or vapor enthalpy, 

is the wall heat transfer for each phase. These terms are left undefined at this level 
because they depend on the specific test case. For each of the test case of section |5j their 
precise expression will be given. 

To analyze what occurs when the volume fraction of one of the phases becomes small, 
this model is too complex. Therefore, in the analysis section below, we focus on the 
isentropic model in one space dimension. 

2.2 Asymptotic analysis of the isentropic two-phase model 

We investigate the behavior of the isentropic two-phase model when the volume fraction 
of one of the phases vanishes. For the sake of simplicity, we consider the one-dimensional 
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model and exclude any source or viscous terms except for the interfacial pressure default 
which makes the system hyperbolic. The isentropic two-fluid model is: 



= 0, 

dtia^PvUv) + dx{a^Pvul) + a^d^p + 5 a^aipul d^a^ = , 
dt{aepeue) + d^^^aipnu]) + atd^p + 5 a^at p ul Dl d^ae = 



Pv 



Pv{p) , 
Pi{p) , 



Civ + ae 



(2.11) 
(2.12) 
(2.13) 
(2.14) 
(2.15) 
(2.16) 
(2.17) 



This model is hyperbolic provided 6 >1 [36] . 

Let us now focus on the behaviour of the system when a phase disappears. We consider 
for instance that the vapor phase is disappearing. The vapor volume fraction becomes 
close to zero. Therefore, it is legitimate to introduce a small parameter e -C 1 which 
measures the order of magnitude of and to rescale follows: 



After this rescaling, the system becomes: 

dt{a^Pv) + d.j,{a^pyU^) = , 
dtiaipe) + d^{aiPiUe) = , 



dt{a^pvUv) + dxia,aPvU^) + dvdxP + ea^atpu^ 5 d^a^ = 
dt{aipeue) + d^^aipeuj) + aAp + e a^aip ul 5 d^ag = , 



Pv 



Pv{p) 
Pe{p) , 



ea^ + = 1 . 

We now write the system obeyed by the formal limit e — )■ 0: 

0, 



dtp I + dx{piUi) = 0, 
dt{a^p.uU^) + dx{a.,pyul) + a^d^p 
dtipeUi) + d^ipiuj) + dxp = 0, 
pv = pv{p) , 
Pe = Pe{p) ■ 







(2.18) 



(2.19) 
(2.20) 
(2.21) 
(2.22) 
(2.23) 
(2.24) 
(2.25) 



(2.26) 
(2.27) 
(2.28) 
(2.29) 
(2.30) 
(2.31) 



Let us make a few comments on the structure of the limit system. First, we notice that 



the system composed of eqs. (2.27), (2.29) and (2.31) is nothing but the isentropic Euler 
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system for a single fluid consisting of the liquid phase. Indeed, the isentropic pressure of 
this fluid is given by the inverse function p{pe) of Pi{p)- The system for the liquid phase 
is thus completely decoupled from the vapor phase. 

Let us now turn towards the system consisting of eqs. (2.26), (2.28), (2.30) which 



determines the vapor variables. Since the pressure p is entirely determined by the liquid 
phase, the pressure term a^dxP in (2.28) is a zero-th order term in a„, multiplied by a 



known coefficient d^p- Therefore, the system for the liquid variables can be written 



So, 



(2.32) 
(2.33) 



where contains only zero-th order terms. The hyperbolicity of the model is determined 



by the left-hand sides of (2.32), (2.33). The corresponding system is a pressureless gas 
dynamics system for the variable U = {avPy,a^pyUy). The pressureless gas dynamics 
system is not hyperbolic. If we write this system dtU + dxf{U) = S, with f{U) = 

(0, Sv), the Jacobian matrix &j does not have a complete basis 



{a^PyUy, ayP^ul) 311(1 S 



of eigenvectors. More precisely, Uy is an eigenvalue of multiplicity 2 but the associated 
eigenspace is of dimension 1. The matrix ^ can be written in the form of a Jordan 
block of size 2. We refer to O [Z| for a detailed analysis of the pressureless gas dynamics 
equations. 

Now, we consider the scaled system (2.19), (2.25). It is a strictly hyperbolic 4 x 
4 system [36j. Consequently, it has a complete basis of 4 eigenvectors. In the limit 
£ — >■ 0, two of these eigenvectors converge towards corresponding eigenvectors of the 
isentropic Euler system for the liquid phase. The other two eigenvectors become parallel 
to each other and parallel to the unique eigenvector of the vapor phase pressureless gas 
system. Appendix A. 2 confirms this deduction: using the first-order approximation of 
the eigenvectors given in [31] for the perfect gas equation-of-state, we show that the 
eigenvectors corresponding to the void fraction and pressure waves become parallel to 
each other when — )■ 0. 

To summarize, this analysis shows that, when a phase disappears, some eigenvectors 
collapse and become parallel. We will see that this phenomenon can raise some issues for 
the numerical scheme. 



2.3 Roe scheme and phase appearance / disappearance 

Roe's approximate Riemann solver [3U [33] is one of the most powerful and widely used 
schemes to solve hyperbolic systems of conservation laws. However, the two-fiuid model 
has non-conservative terms. Toumi and Kumbaro [36] have proposed a generalization 
of the Roe linearization to non-conservative systems. The non-conservative two-phase 
system can be written in the quasi-linear form: 

-+A(V)-^0. (2.34) 
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In the finite volume framework, tlie generalized Roe scheme can be written as: 



yn+l _ yn 1 

' ' + ^(<f-(V„V,+i) + '^'^(V.-i,V,)) =0, (2.35) 

with 

$± ( Vi, V,+i) = A± ( 1 ) ( V,+i - Vi) . (2.36) 

The Roe matrix is the Jacobian matrix A of the system taken in an appropriate lineariza- 
tion state Vj^i. For a non-conservative system, the linearization state is chosen so that 
the shock waves at the interface between cells i and i + 1 remain those of an equivalent 
conservative system [SB!- The positive and negative Roe matrices are defined by: 

A(V,, i) ± |A(V., i)| 
A^(V,^.) = ^ . (2.37) 

where |A(V-_^i)| is the absolute value of A(V^_,_i). In the OVAP code, the second term 
in (2.35) is evaluated implicitly and the resulting nonlinear system for Vl^^^ is solved by 



Newton's iterations 

The computation of the absolute value of |A(Vj_,_i)| is performed as follows. Let A be 
a diagonalizable matrix. We write 

A = 7^diag(Al,...,A7v)7^"^ (2.38) 

where the 's are the eigenvalues of A, diag(Ai, . . . , Aat) is the diagonal matrix whose 
diagonal coefficients are the A^ 's, and TZ is the matrix whose columns are the eigenvectors 
of A. Then, |A| is given by 

|A| =7^diag(|Al|,...,|A^^|)7^-^ (2.39) 



Formula (2.39) for the matrix absolute value is valid as long as A is diagonalizable. 
However, if the system loses its hyperbolicity, the eigenvectors of the Jacobian matrix A do 
not form a complete basis anymore, the matrix TZ becomes singular and strictly speaking. 



|A| is no more defined. We have seen that the limit system (2.11 )-(2.17) is not hyperbolic 
for a„ — )■ for the precise reason that the eigenvectors do not form a complete basis any 
longer. In practice, during a computation, numerical problems begin to appear with the 
Roe scheme for a„ G [10~^, 10""^], depending on the considered case. These problems are 
caused by some of the eigenvectors becoming almost parallel when the volume fraction 
decreases. The matrix TZ becomes highly ill-conditioned. The numerical accuracy of the 
eigenvector decomposition is then strongly affected. Therefore, the use of the Roe scheme 
based on an eigenvector decomposition of the Roe matrix must be avoided when phases 
appear or disappear. We will see that |A| can be computed with different methods which 
do not require the use of the eigenvector decomposition of A. With this aim, we recall a 
certain number of results stemming from functional calculus 



8 



Let A be a diagonalizable matrix and denote by Sp(A) = {Ai, . . . , Ajv} the spectrum 
of A. Let $ be a continuous function defined on an open interval X containing Sp(A). 
The matrix $(A) is defined by 

$(A) =7^diag($(Al),■■■ ,$(Aiv))7^-^ 



with TZ defined by (2.38). We note that $(A) only depends on the values of $ on 
Sp(A). Additionally, if is a sequence of function such that ($„(Ai), . . . , $„(AAr)) — )■ 
($(Ai),. . . ,$(A7v)) in M^, then $„(A) — )■ $(A) in any matrix norm. Of course, this is 
the case if ||<l>„ — $||oo — 0. Here, ||$||oo denotes the uniform norm in the space C^{X) of 
continuous functions on the closure X of I. 

Consequently, if a function approximates the absolute function the resulting 
$(A) approximates |A| to the same order. Thus, we are looking for approximation func- 
tions $ which allow the computation of <I>(A) without requiring the eigenvector decompo- 



sition (2.38). This can be achieved by taking $ as a polynomial P such that -P(Aj) ~ |Aj|, 
for all z = 1, . . . , iV. Indeed, -P(A) can be simply calculated by taking successive pow- 
ers A'^ of A and does not require the eigenvector decomposition. This gives rise to the 
so-called polynomial schemes [10]. Then, we will also consider an alternative, consisting 
in using the hyperbolic tangent function, which can be computed by solving a matrix 
ordinary differential equation. In all these cases, $(A) will still be defined even when A 
ceases to be diagonalizable and the scheme will not breakdown at phase appearance or 
disappearance. 

2.4 Eigenvalues of the full two-phase model 

Although the full two-phase model of section |2.1| is not as simple as the isentropic model 



of section 2^, some information about the eigenvalues and eigenvectors of the system 
can be obtained. Because of the complexity of this model, no analytical expression of the 
eigenvalues is available. However, approximations given in [251 ES] enable us to discuss the 
behavior of the eigenvalues when a phase appears or disappears. The detailed computation 



is given in appendix A.l 



Since the hyperbolicity of the model is only determined by the left-hand sides of eqs 



(2.1)- (2.6), the precise knowledge of the source term is again unnecessary. The system is 
posed in dimension d. So, there are 4 + 2(1 eigenvalues of the Jacobian matrix. In general, 
there are two fast eigenvalues which are of the order of ± c, where c is a characteristic 
sound velocity of the two-phase mixture, two eigenvalues of the order of called the void 
eigenvalues, and two trivial eigenvalues, each of multiplicity d, respectively equal to the 
vapor and liquid velocities and U£. Note that the void eigenvalues can be complex if 
the interfacial closure terms are not carefully chosen (see [28]). The fastest eigenvalues 
Uy ± c are always real and remain distinct from the other eigenvalues. We will denote 
them by Xmax for the largest and A^m for the smallest. All the other eigenvalues, that 
we will call the " intermediate eigenvalues" and collectively denote by ;\j5'termediate j^g^yg i^j^g 
same orders of magnitude as long as the two-phase flow stays subsonic. In the example of 
a boiling channel which is relevant for our applications (see section |5]), the ratios between 
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the orders of magnitude of the fastest eigenvalues and the intermediate eigenvalues are 
the following: 



I \ intermediate I I \ intermediate I 



— — - ^ 10"^ (2.40) 

I ^max I I '^min \ 

Suppose now that the vapor phase disappears : a^, — j- 0. Then, the fast eigenvalues 
tend towards ± where the expression of is given in appendix A.l They remain 
distinct and the associated eigenvectors do not collapse. However, the void eigenvalues 
become of the order of magnitude of and the corresponding eigenvectors collapse. 
Qualitatively, the same phenomenon as in the isentropic case occurs: two eigenvectors 
become parallel in the limit a^, — )■ and the eigenvectors do not form a complete basis 
any longer. The matrix TZ formed by the eigenvectors becomes ill-conditioned. The 
computation of the Roe matrix becomes highly inaccurate. 



3 Numerical schemes based on polynomial or hyper- 
bolic tangent evaluations of A 

As already announced, polynomial schemes avoid the use of the eigenvector decomposition 
of the Roe matrix A to compute |A|. In this section, we give a presentation of polynomial 
schemes and provide a selection of high-degree polynomials which are well-suited to multi- 
phase flow calculations in the situation of phase appearance or disappearance. Polynomial 
schemes have been introduced in [10] and used in [2H [28] . We also present an alternative, 
based on the evaluation of A using the hyperbolic tangent function. This method is, to 
the best of our knowledge, new. 



3.1 Computation of |A| with a polynomial 

We recall the approach sketched at the end of section |2.3[ It relies on the approximation 
of |A| by a polynomial P such that -P(A) ^ |A|. Indeed, the matrix polynomial 

n 

P(A) = ^a,A^ (3.1) 

k=0 

of a diagonalizable matrix A can be alternately computed, using the eigenvector decom- 
position (2.38), by: 

P(A) = 7^diag(P(Al), . . . , P(A^)) n-\ (3.2) 

Therefore, if P satisfies 

P(A,) = |A,|, VzG[l...Ar], (3.3) 
i.e. if it interpolates the absolute value function at all the eigenvalues of A, then 

P(A) = |A|. (3.4) 
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Therefore, there are two ways of computing |A|: either by formula (2.39), or by (3.1) with 



a polynomial P satisfying (3.3). However, the advantage of formula (3.1) over (2.39) is 



that it does not use the eigenvector matrix TZ and is consequently faster. Additionally, 
the computation of P{A) does not breakdown if the matrix is not diagonalizable, while 
that of |A| does. In fact, |A| is no more defined in this case while P{A) stays defined. 
Therefore, the polynomial formula for |A| is better suited to the case where the eigenvector 
decomposition of A breaks down and the matrix TZ becomes ill-conditioned. In view of the 



discussion of section 2^, polynomial schemes appear as methods of choice for situations 
of phase appearance or disappearance 

In practice, the selection of the polynomial is crucial. Indeed, it may be useless to 



verify (3.3) exactly, i.e. for all the eigenvalues. It may increase computational costs to 



no avail and may be detrimental to the stability of the scheme. If (3.3) is not satisfied 
exactly, then P(A) f» |A| instead of satisfying (3.4) exactly. The selection of P becomes a 



compromise between accuracy on the one hand, and stability and computational efficiency 
on the other hand. We discuss these issues below. 

For explicit schemes, Degond and al [lOJ have shown a sufficient stability condition 
for polynomial schemes, under the CFL condition. Let Amin and Amax be the smallest 
and largest eigenvalues of the Roe matrix A, and Omax = niax{|Amin|, |Amax|}- Then, the 
stability criterion reads 



\x\ < P(x) < an 



Vx G [Amm; Aj^ 



(3.5) 



This condition is represented graphically in FigjTJ The graph of the polynomial in the 
interval [Xmin, ^max] must be contained in the coloured area of the figure. 



Figure 1: Stability condition. The graph of the polynomial in the interval [Xmin, Xmax] 
must be contained in the coloured area in order to ensure the stability of the explicit 
scheme. 



Condition (3.5) ensures the stability of the scheme. Accuracy requires that large 
oscillations of the polynomials near the eigenvalues should be avoided. Indeed, if the 
derivative of the polynomial about one of the eigenvalues is large, round-off errors may 
be amplified. A small difference between the true eigenvalue A and the computed one A 
may cause a huge discrepancy between |A| and -P(A), thus creating numerical inaccuracies. 
In [2U [2S], the Lagrange interpolation of the |Aj| 's is used in the Newton basis. This 



interpolating polynomial, further on referred to as 'Pexact' verifies (3.3) exactly but has 



large oscillations. The resulting method is as good as the classical Roe scheme in standard 
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(a) Poxact (b) Zoom on small eigenvalues 

Figure 2: Exact interpolating polynomial Pexact- The black line is the polynomial, the red 
line is the absolute value function, and the blue spots are the eigenvalues. Left, the graph 
of the polynomial in the full range of interest [Xmin, ^max]- Right: a blow-up of the graph 
in the region of eigenvalues of the smallest absolute values. 



situations. However, it breaks down at phase appearance or disappearance due to the large 
oscillations that are generated at the extremal eigenvalues (see fig. |2]). These oscillations 
are caused by the intermediate eigenvalues which get very close one to each other (see 
section [2]). This induces ill-conditioning of the matrix involved in Newton's method of 
computation of the Lagrange interpolation polynomial and loss of accuracy. In what 
follows, we develop new approximating polynomial avoiding this difficulty. 

In [lOJ, a method is presented to approximate |A| using interpolating polynomials P„ 
of degree n = 0, 1 or 2 (resp. denoted by Pq, Pi, P2). The interpolation only focuses 
on the extremal eigenvalues Xmin and Xmax, and adds a condition over one derivative 
in the P2 case. Fig. |3] depicts the graphs of the Pq, Pi and P2 polynomials. They all 



respect the stability condition (3.5). For the targeted applications in which the orders 



of magnitude of the eigenvalues satisfy (2.40), it appears that the absolute values of the 



intermediate eigenvalues, which are close to zero, are not approximated accurately enough. 
This results in a quite poor approximation of |A| and gives rise to very diffusive schemes. 
These schemes are not accurate enough and will be discarded. In the following sections, 
we propose the construction of new polynomials which considerably improve the accuracy 
while maintaining the stability of the scheme. 

3.2 Approximation of |A| by high-degree interpolating polyno- 
mials 

We look for polynomials that provide accurate approximations of the absolute value func- 
tion on the spectrum of the matrix, while maintaining the stability of the scheme when a 



phase appears or disappears. Such polynomials must satisfy the stability condition (3.5) 



and avoid large oscillations near the eigenvalues. To meet the accuracy constraint, we need 
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(a) Po 



(b)Pi 



(C) P2 



Figure 3: Interpolating polynomials Pq Pi and P2 based on the extremal eigenvalues only. 
They respect the stability condition (the stability area is coloured). 



to consider high-degree polynomials. We consider polynomials interpolating the absolute 
value function in the interval [—1, 1]. The general case can be deduced by applying the 
result to the matrix A/a^^^, with a^ax = max^ |Afc|. Two approaches are considered: fixed 
interpolation and dynamic interpolation. Fixed interpolation means that the approximat- 
ing polynomial does not depend on the eigenvalues and approximates the absolute value 
function in the whole range [—1,1]. Dynamic interpolation means that the approximating 
polynomial depends on the eigenvalues to be interpolated and focuses on the quality of 
the approximation near these eigenvalues. The second approach, although slightly more 
time-consuming, since it requires to re-construct the polynomial at each time-step and 
each cell- interface, will reveal to be more efficient. 



3.2.1 Fixed polynomial interpolation 

Polynomials interpolating extremal points: the P2p polynomials. A first idea is 
to construct even polynomials P2p = J2k=o '^kx'^'' p E N. P2p is constructed such that 
P2p{'x) vanishes at a; = 1 as well as its derivatives up to the order p. The polynomial 
is then uniquely defined by: 



P2p is even. 



^2p(l) 
^2.(1) 



1, 



(3.6) 



3 



,P- 



The coefficients ak are calculated once for all by solving a linear system. The larger the 
order of contact of P2p{x) — \x\ with 0, the better the approximation is. Fig. 4(a) displays 



P2p{x) for p = 1 to 16. The value p = 1 corresponds to the interpolating polynomial 
P2. Higher values of p clearly provide better approximations of However, due to 
the inversion of the linear system, calculating the coefficients ak beyond p = 16 presents 
numerical instabilities. The so-obtained accuracy is still not entirely satisfactory. 



Even polynomials Prdf interpolating intermediate points. In this example, the 
even polynomial P{x) = Yl^=o^kx'^'^ of degree 2p interpolates |a;| at a series of m points 
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(a) Even polynomials P2p. 



(b) Phdf 



Figure 4: (a) Even polynomials interpolating the absolute value at the extremal 
points of the interval [—1, 1], for p = 1, . . . , 16. (b) Even polynomials Phdf interpolating 
intermediate points. 



{xi)i<i<m, Xi e [0,1], such that P{x) — \x\ has a contact of order q with at Xj. The 
polynomial is determined by: 



P2p is even, 

P{Xi) \Xi\ 
P\X,) = 1 

pO)(a;,) = 0, J =2 



(3.7) 



, . . . , 



Ci 




The degree of the polynomial is 2^™j^(q + 1). After several trials, it appeared that an 
almost optimal choice was obtained with the following parameters: 



(3.8) 



where the ti are the two Tchebychev points on the interval [0, 1]. The Tchebychev points 
give minimal oscillations for an interpolating polynomial, and have the following expres- 
sion for an interval [a, h] divided in n points: 

h — a (2k + 1)71 a + b ,„ 

tk = cos^ V + 3.9) 

2 2(n + l) 2 ^ ' 

We will refer to this polynomial as Phdf, for "High-Degree Fixed" polynomial. Its coeffi- 
cients are calculated once for all by solving the linear system (3.7). The numerical values 
of the coefficients are given in appendix A. 4 As we have 

min {Phdf{x) - \x\) ~ -lO-^^^ 

x&[0,l] 
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we add a constant equal to 10 so that the polynomial remains greater than the absolute 
value function. The so-obtained polynomial respects the stability condition. The graph 



of Phdf is given in Fig. 4(a) (b). The approximation of the absolute value function by 
Phdf is improved. However, the absolute value of the intermediate eigenvalues, those 
which have a magnitude close to 0, remain inaccurately approximated and the scheme 
remains too diffusive. 

3.2.2 Dynamic polynomial interpolation 

As seen is the previous section, to improve accuracy, it is necessary to take into account 
the intermediate eigenvalues. The resulting polynomial depends on the eigenvalues, and 
motivates the terminology 'dynamic interpolation'. 

One of the reasons for the large oscillations of the interpolation polynomial Pexact 
is the presence of a cluster of intermediate eigenvalues near 0, which are very close to 
each other and which are 4 orders of magnitude smaller than the extremal eigenvalues 



(see (2.40)). In [23], an approximate Roe matrix is generated by treating this cluster of 
eigenvalues like a single eigenvalue (the largest of them). The method of |23j provides 
comparable results to the usual Roe method in standard situations. Inspired by this work, 
we generate an approximate polynomial which interpolates the extremal eigenvalues \mini 
Xmax and only one of the intermediate eigenvalues, the largest one, denoted by A^^^ (as 
well as its opposite — A^"^ for symmetry reasons). It thus avoids the interpolation of 
the cluster of very close intermediate eigenvalues. Conditions on derivatives are added so 
that the stability condition is respected locally around the eigenvalues. The polynomial 
is a Hermite interpolation polynomial constructed in Newton's basis, and is calculated 
at each time-step and for each cell-interface. The computation does not break down at 
phase appearance or disappearance. Indeed, the collapse of eigenvalues only concerns 
intermediate ones. The design of the polynomial considers already only one of them and 
it does not matter how many of them are distinct. 

We will call this polynomial Phdd, for "High-Degree Dynamic" polynomial. It verifies 
the following conditions : 



max I ; 



I, PHDoi^max) = |A 

^HDD[±\nt ) — Wnt 1 5 ^HDD\^min) — — J-, 

PHDoi^max) = 1, PRDoi^'^inT) = "^5 

PHDoi.'^lnt^) ~ 1? PHDD^^min) = PuDD^^rnax) =0, j = 2, 10. 



As fig. 5(a) (left) shows, the contact between the polynomial and the absolute value is 
very good in the neighborhood of the extremal eigenvalues, as the first derivative is set 
to ±1 (the derivative of the absolute value) and the other derivatives are set to zero. The 
large oscillation between the extremal eigenvalues and the intermediate eigenvalues is not 
a problem because there exist no eigenvalues in this region. This allows us to approximate 



x\ near x = very accurately. Fig. 5(b) (right) shows a blow-up close to zero. We notice 



that the stability condition (3.5) is indeed verified for all intermediate eigenvalues, as we 



do have Phdd{,X) > |A| for all intermediate eigenvalues. 
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(a) XZT = 10"^ (b) Zoom on x = 

Figure 5: Dynamic interpolating polynomial Phdd based on the interpolation of the ex- 
tremal eigenvalues and the largest intermediate eigenvalue (left). There are no eigenvalues 
in the region of the large oscillations in the median regions and the stability condition 
(3.5) is respected locally about the eigenvalues. The right figure shows a blow-up of the 
left picture near 



We will test the behavior of the Ph d f and Ph d d polynomial schemes in the numerical 
section |5| We will see that numerical difficulties are reduced but some positivity problems 
remain. In the following part, we develop a method to specifically treat the positivity 
problems. It will rely, among others, on the possibility of tuning the amount of diffusion 
in the Phdd polynomial. But first, let us detail an other way to compute |A| without 
using the eigenstructure of the matrix. Based on the same principle as the polynomial 
solvers, the method uses the hyperbolic tangent function. 



3.3 Approximation of |A| by means of the hyperbolic tangent 

In this section, we present an alternative to the use of the polynomial schemes. We recall 
that the goal is to compute the absolute value matrix |A| without using the eigenvector 
decomposition of A. We introduce the following approximation $(x) of the absolute value 
function 

$(a;) = r + (1 - r) X tanh(-) cotanh(-). (3.10) 

r r 



with 



tanh(x) = ^ , cotanh(x) 



+ tanh(a;) 
and r > is a parameter. 



As in section |3.2[ we will normalize the matrix A by the largest absolute value of the 
eigenvalues and study the function $ only in the interval [—1,1]. We have $(1) = 1 and 
$(x) > |a;|,forallx G [—1,1]. Furthermore, it is easy to realize that | tanh(^)— S'ign(x)| — > 
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0, and consequently, that $(a;) — I — > when r — )■ 0, uniformly for x G [—1, 1], where 



S'ign(x) denotes the sign function. Therefore, $ obeys the stability condition (3.5) (see 
also fig. [T| and is an approximation of with a controlled accuracy. The graph of ^ is 
represented on fig. [6] for different values of the parameter r. 
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(a) T = 0.1 



(b) T = 0.001 



Figure 6: The function $(x) (in red) as a function of x G [—1, 1]: the parameter r controls 
the accuracy of the approximation of |x| (in blue) by $(x). 

The practical choice of r is performed as follows. Our goal is to approximate closely 
all the non-zero eigenvalues of A including the smallest. Consequently, we must have 
T < mm|Aj| = Xs- Tab. jlj presents the maximum error on the interval [A^, 1] between 

|x| and $(x) for different values of r and for = 10"^. The choice of r = A^/IO gives 
already quite good accuracy. Nevertheless, if we want to add some diffusion, it is possible 
to reduce the accuracy by increasing the parameter r. 



r 


max $(x) — X 


10 " 
100 ^ 

= 10"'^ 

1000 


9.998 X 10-*^ 

9.998 X 10-^ 

9.999 X IQ-^ 



Table 1: Accuracy of the approximation of the absolute value function depending on r 

We now present how we compute the hyperbolic tangent of the matrix A without 
using the eigenvectors of the matrix. We note that the scalar hyperbolic tangent function 
X — )■ tanh(ax) (where a G M is a constant) satisfies the following differential equation: 



d 

dx 



tanh(Q;x) = a{l — tanh(ax) 



(3.11) 



Therefore, we solve the matrix differential equation: : 

rfx(C) 



dc 

X(0) 



(3.12) 



0. 
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which yields = tanh(^A) with ^ G M. We solve this differential equation by means 
of an iterative implicit method. An explicit method based on a fourth order Runge-Kutta 
method can also be used but the computational cost is higher. The scheme is written: 



where h is the iteration step. We use the Newton method to find X'^^^ at each iteration. 
The number of steps needed to prevent the algorithm from diverging can be set to a 
constant value Ni. Each iteration contains another loop of maximum N2 iterations to 
find X*^"*"^ by the Newton method. We have chosen A^i = 100 and = 40. 

4 Numerical treatment of positivity losses 

Positivity problems tend to appear during the simulation of phase transitions. We will 
first review previously developed positive schemes. We will then propose an other method 
to solve the positivity problems, method complying with the constraints mentioned in the 
introduction : no analytical expression of the eigen-elements are available, the computa- 
tion of the eigenvectors should not be used during phase transitions, and the method has 
to be compatible with an implicit scheme or large time-steps. We will give the general 
principle of the method and then the features of its implementation. 

4.1 Previous works on positive numerical schemes 

A positive scheme for the two-fluid two-phase flow model has been proposed in [8]. The 
explicit scheme introduces a splitting in the resolution of the bifiuid model. The first step 
(hydrodynamic step) solves separately two uncoupled full Euler systems, for each phase, 
by means of a kinetic solver, for stability reasons. The non conservative terms in d^a are 
reformulated and included in the source terms. A second step enforces the equality of the 
pressures and allows to compute directly the void fraction and the pressure. This scheme 
preserves the positivity of all thermodynamics variables under a CFL-like condition. 

As few works concern the resolution of the two-fluid two-phase flow model specifically, 
let us also mention the positive numerical schemes that have been designed for Euler 
equations, or gas dynamics equations. 

Einfeldt and al. [H] consider the HLLE solver for the Euler equations. HLLE is 
positively conservative, but less accurate than the Roe scheme. Anti-diffusion parameters 
are introduced in the HLLEM scheme to take out excessive dissipation. In |30] , Perthame 
and Shu show that the Lax-Friedrichs scheme is positivity preserving, which echoes the 
general observation that the more diffusive a scheme is, the more robust it is, robustness 
including here positivity preservation. 

For the sake of accuracy, other works introduce second order schemes. The positivity 
constraint is ensured by limiting techniques. In [21j , the symmetric limited positive scheme 



x'^ ^ x(a) 

X(0) = 



(3.13) 
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(SLIP) conserves the positivity thanks to the use of a hmited diffusive flux which makes 
the scheme local extremum diminishing (LED). This property is stronger than the total 
variation diminishing (TVD) property proposed by Harten [18] (if the scheme is LED, then 
it is TVD, LED and TVD being equivalent in one dimension), and ensure that a local 
maximum cannot increase, and a local minimum cannot decrease. Thus, if the solution 
is positive at one moment, than the global minimum is positive and cannot decrease and 
become negative. This SLIP scheme can be applied to the Roe scheme for a system 
of conservation laws. The construction of the scheme requires the computation of the 
eigenvectors. 

The Maximum Limited Gradient reconstruction technique [2], also gives a second 
order positivity preserving method. In [27], Liu and Lax propose a family of second order 
positive schemes for mult i- dimensional hyperbolic systems of conservation laws, using a 
limiter in the numerical flux. The eigenvectors and eigenvalues of the Roe matrix have 
to be evaluated explicitly to construct the scheme. The second order central scheme 
described in [22] is based on a predictor corrector method and is positive due to the scalar 
maximum principle. An other positive scheme based on flux limiters can be found in |3]. 

Other methods designed to address the lack of positivity have also been based on a 
modification of the Roe scheme. Dubroca in [12] and Gallice in [16] propose extensions 
of the Roe's solver for the Euler equations on the one hand and for systems of magne- 
tohydrodynamics on the other hand. In [TJ], Einfeldt and al. concluded that there is 
no positively conservative Roe matrix, but also specified that this statement only applies 
to Roe's matrices based on Jacobian matrices. These works introduce the derivative of 
the pressure in the direction of the fluid velocity in the Roe matrix. This decomposition 
allows to introduce parameters that can be chosen so that the solver becomes positively 
conservative. The demonstrations of the positivity of the two methods are based on the 
eigenvalues and eigenvectors of the Roe matrix, that can be calculated analytically for 
the considered systems. 

An noticeable point in [TB] is the establishment of a link between diffusion and pos- 
itivity. Gallice introduces parameters in the Roe matrix as a way to exactly tune the 
dissipation so that the scheme becomes positive. Our work is based on the same idea: 
finding the right amount of diffusion so that the positivity problem can be solved. 

Up to now, most of the methods available in the literature are based on the construc- 
tion of specific schemes which are designed to prevent the positivity problems. A different, 
less developed strategy is to treat the positivity problem when it appears rather than try- 
ing to prevent it. Such a strategy is proposed by Romate in [32]. Romate remarks that 
the HLL scheme [13] is positive under a CFL condition but too diffusive for practical use 
as such, while the Roe scheme is accurate but is not positive. He presents a combination 
of the two schemes. The Roe flux is used until the appearance of a positivity loss in an 
adjacent cell. If such a problem occurs, the time-step computation is restarted using the 
the HLL flux. The newly computed cell value will therefore be positive. 

The method presented below is inspired from Romate's strategy [32]. Rather than 
trying to prevent positivity problems from appearing, it consists in applying a special 
treatment when such problems appear. It is also inspired from Gallice [TB] in that the 



19 



treatment consists in increasing the numerical diffusion in some way. 



4.2 Description of the method 

The method is inspired from [12] where positive Roe schemes for the gas dynamics and 
MHD equations are proposed. This work is based on the fine tuning of the numerical 
diffusion in order to make the scheme positive. The availability of analytic expressions for 
the eigenvalues of the Roe matrix is a key ingredient in the demonstration of the positivity 
preserving property. In our case, the two-phase model is too complex to allow for analytic 
expressions of the eigenvalues. Thus, we will develop the idea in a different way. We 
increase the numerical diffusion where positivity losses are detected. To this purpose, 
we use the polynomial scheme based on the Phdd polynomial just described. Indeed, 
it provides an easy way to adjust the numerical diffusion. We stress that, by contrast 
to [16], we have no proof that positivity is preserved, but only numerical evidence that 
robustness against loss of positivity is enhanced. 

To introduce numerical diffusion within the Phdd polynomial, we just reduce the accu- 
racy of the interpolation of the intermediate eigenvalues, thus making the scheme more dif- 
fusive. Instead of the points (itA™"^, |A^j^|), we interpolate the points (±A^"^, \D A™j^|), 
where D > 1 is a diffusion coefficient which is chosen to maintain the condition P(A™"^) > 



|A™"^|, in agreement with the stability condition (3.5). 

When no positivity problem appears, the code is normally run with the Phdd polyno- 
mial associated to a diffusion coefficient D = 1. If, after a time-step, a positivity problem 
is detected in some cell, the computation of the time-step is restarted using a new Phdd 
polynomial using D = 10 in the adjacent cell interfaces. If the positivity problem remains, 
we further increase D (the precise algorithm is given below) . If finally, the maximal value 
of D is reached (beyond which D\X^^^\ > ma.x{\ Xmin\, ^max), which is forbidden by the 
stability condition ( |3.5 )), then, we stay with this value of D and reduce the time step. 
The details of the algorithm are given below. 

In the numerical section [5} we will see that the combination of the Phdd polynomial 
scheme and the present treatment of positivity problems leads to significant improvements. 
The robustness and reliability of two-phase flow codes in situations of phase appearance 
and disappearance is greatly enhanced. 



4.3 Implementation 

We use the Phdd polynomial scheme to adapt the diffusion by the means of the coefficient 
D so that we interpolate D |A™"^| instead of |A^"^|, the largest intermediate eigenvalue. If 
D > 1, the diffusion is increased. A coefficient Cj is attributed to each cell to take inventory 
of the occurrence of positivity problems within a given time-step. The treatment proceeds 
according to the following algorithm, which describes a time-step advance. We will call 
this method Ph°dd- 

a - At the beginning of the time-step, the counter q = on all cells: no positivity 
problems has occurred yet. 



20 



b - On all interfaces, compute \A\ with Phdd and D = 1. 
c - Solve the time-step. 

d - Loop around all the cells to check for problems. For all cells i where a positivity 
problem or a convergence problem during the computation of the variables of states 
(pressure, enthalpy) has occurred, increase the counter q by 1. 

e - Compute \A\ again with an increased diffusion on the interfaces where at least one of 
the neighboring cell is such that q 7^ 0. In order to increase progressively the diffusion, 
we use D = lOcf. IfD |A^?^| > max(|A™„|, A^^x), we set D = max(|A^,„|, A„ax)/|A™r' 



in order to remain in the domain of validity of the stability condition (3.5). 



f - Solve the time-step again and iterate if necessary until all problems are solved. We 
set up initially a maximal number of iterations. If this number is reached, we reduce 
the time-step At and restart the computation of the time-step. 

This method allows to overcome most positivity issues, but also problems which may 
occur in the computation of the variables of states (pressure, enthalpy). It is very im- 
portant to note that, as diffusivity is added very locally, this method has no negative 
repercussion in terms of global accuracy, which is preserved. Let us also note that, in 
|16j . there is no analytic expression of the parameters introduced to correct the system 
matrix so that the scheme is positive. These parameters have to be "large enough" to 
ensure the positivity, and their value is also obtained by an iterative procedure, as in our 
method. 



5 Numerical results 

We present several test-cases in one and two dimensions. The Ransom faucet test-case will 
first allow us to compare the accuracy of the different polynomial schemes. The boiling 
channel in the saturated case will then highlight the improvement brought by polynomial 
schemes in a situation of phase appearance or disappearance. We will then show more 
difficult test-cases where the positivity treatment is required : the boiling channel in 
the subcooled case, and the two-dimensional tee-junction test-case. In all test-cases, the 
water-and-steam equation of state will be used. The International Association for the 
Properties of Water and Steam (lAPWS) provides internationally accepted formulations 
for the properties of steam and water. In all test cases, the used model is that of section 



2.1, and the right-hand sides will be specified precisely. 



5.1 Ransom faucet 

This non-stationary one-dimensional test-case was proposed by Ransom in [T3]. It con- 
siders the flow of a water column at the outlet of a faucet opening out into a vertical 
enclosure containing air. The considered tube is 12 m high, and the inlet velocity is 10 
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m/s whereas the air is at rest. The ratio of the sections of the nozzle and of the en- 
closure is such that the integrated void fraction over the section is equal to 0.2. In this 
configuration, a striction phenomenon of the jet is observed due to the effect of gravity. 
Indeed, if we make the assumption that the jet remains coherent (no tear-off of liquid in 
the form of drops, no penetration of air into jet), the acceleration of the liquid due to 
gravity necessarily results in a narrowing of the cross section of passage of the liquid, by 
conservation of the fiowrate. Moreover, as the initial conditions correspond to the solution 
which would be obtained in the absence of gravity (therefore = 0.2 everywhere), a void 
fraction discontinuity is propagated from the inlet section to the outlet section as from 
the initial time. 

This test-case allows to evaluate the accuracy of a scheme, the amount of numerical 
diffusion being visible on the void fraction front. The velocity of the void fraction wave 
has to be correctly captured. The model is the six equations two-fluid model presented in 
section 2.1, The only source terms are the interfacial pressure term (2.10) with 6o = 1.1, 
and the gravity with g = 9.81m/s^. At the inlet, the following values are fixed: = 0.0 
m/s, Ui = 10.0 m/s, K = 324.594 kJ/kg, hg = 209.283 kJ/kg, and = 0.2. At the 
outlet, the pressure is fixed: Poutiet = 10^ Pa. The computational method is explicit. We 
used a one- dimensional mesh with 100 cells. The maximum time is t = 0.6 s. 

Fig. [7| represents the profiles of the volume fraction, pressure and velocities, for the 



Phdf and Phdd polynomial solvers, and the Tank scheme described in section 3.3 The 



results are compared to the solution given by the Roe scheme, and in the case of the 
void fraction, the analytical solution is shown. We can see that the high-degree dy- 
namic interpolating polynomial Phdd has an equivalent accuracy as the Roe scheme. 
The high-degree fixed interpolating polynomial Phdf shows good accuracy. This was 
to be expected in this test the vapor and liquid velocities are large. Thus, the 

intermediate eigenvalues whose orders of magnitude are the fluid velocities are in a range 
where Phdf approximates accurately the absolute value function, yielding an accurate 
result. The Tank scheme also has an equivalent accuracy as the Roe scheme but, due to 
the computation of the hyperbolic tangent of a matrix, the computational cost is high: 
the computation cost is 145 times larger for the Tank than for the Phdd scheme. 



5.2 Boiling channel 

The test-case consists in a one-dimensional vertical channel of length Lh = 3.65m with 
upward flowing water [H [37j. A uniform heat flux is imposed along the wall of the channel 
and causes the appearance of vapor. Two cases are considered: at the entrance, the water 
can be either saturated in vapor or be subcooled (i.e. be colder than the saturation 
temperature where vapor starts to form). In the first case, vapor creation starts at the 
inlet. In the second case, vapor creation starts further in the channel, when the saturation 
is reached, for y = ybou- This point can be estimated analytically and is yhou = 1.21m, 
with the data used in the present test-case. This test-case checks the ability of the scheme 
to deal with a large range of volume fractions and to capture the onset of boiling yi^ou in the 
subcooled case correctly. The physics includes stiff source terms and couples hydraulics 
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Figure 7: Ransom faucet with 100 cells. Void fraction (a), pressure (b), liquid velocity (c) 
and vapor enthalpy (d) at time t = 0.6 s, as functions of the height in the column. The 
Roe scheme (red squares), Phdf (green diamonds) and Prdd (blue triangles) polynomial 
solvers, and the Tank (violet triangles) method are compared. The analytical solution 
for the void fraction is shown in black dashed line on figure (a). 



with wall heating. 



The model is the six equations two-fluid model presented in section 2.1 Physical 



sources include drag force, wall friction, mass and heat transfer, and gravity. The detailed 
expressions of these source terms are given in appendix A. 3.1 The computation is implicit. 
We used a one-dimensional mesh with 150 cells. We show the results at t = 5s. At the 
inlet, the following values are fixed: Uy = 0.7802 m/s, ui = 0.7802 m/s, hy = 2.784e6 
kJ/kg. The liquid enthalpy is hi = 1262 kJ/kg in the saturated case and hi = 1029 kJ/kg 
in the subcooled case (it corresponds to a subcooling of AT = 45°C, i.e. the temperature 
is lower by 45°C to the saturation temperature at which vapor starts to appear). The inlet 
fluid is supposed to be pure water. Thus, the initial and inlet vapor volume fractions aj, 
will be set as small as possible. At the outlet, the pressure is fixed: Poutiet = 68.73 10^ Pa. 
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Figure 8: Boiling channel in the saturated case with 150 cells. Void fraction for a„ = 10~^ 
(a) and for = 10~* (b) at time t = 5 s as a function of space. Results by the Phdf 
(green squares) and Prdd (blue diamonds) polynomial schemes, and the Tank (violet 
triangles) method. The Roe scheme (black circles) is depicted on fig. (a) but breaks 
down in case (b). 



5.2.1 Boiling channel: saturated case 

In the saturated boiling channel test-case, the heating sparks the creation of vapor from 
the inlet of the channel. The volume fraction range goes from zero to 0.95. In practice, 
we will try to set the initial and inlet vapor volume fraction as close to zero as possi- 
ble. This case is a good demonstration of the relevance of polynomial schemes for phase 
appearance or disappearance. Indeed, the standard Roe scheme breaks down when vapor 
volume fractions are smaller than = 10~^. The polynomial schemes Phdf and Phdd, 
and the Tank scheme, have no problem whatsoever even for volume fractions as small as 
ai = 10-«. 

The profile of the vapor volume fraction is shown on fig. |8} As the flow is saturated, 
the vapor volume fraction starts increasing at the inlet of the channel. We can see on 



fig. 8(a) (a) for = 10~^ that the Phdd and Tank methods have the same accuracy 
as the Roe scheme, while Phdf is significantly more diffusive. This is due to the large 
discrepancy between the extremal and intermediate eigenvalues, which are of the order 



of magnitude given by eq. (2.40). In this case, the intermediate eigenvalues are not 



approximated very accurately by the polynomial Phdf and the result is diffusive. We 



show on fig. 8(b) (b) the void fraction profile obtained by the Phdf, Phdd and Tank 



schemes for a\, = 10 . The positivity treatment is not needed on this case. 



5.2.2 Boiling channel : subcooled case 

In the subcooled boiling channel case, the vapor starts being created when the saturation 
is reached, at the boiling point ybou = 1.21m. This case is more difficult than the saturated 
case because fluctuations are created at the boiling point and positivity problems often 
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occur at this position. On this case, the Roe scheme presents problems for void fractions 
smaller than aj, = 10~^. The profile of the vapor volume fraction for an inlet and initial 
vapor volume fraction of = 10~^ is displayed on fig. 9(a) for the Roe scheme and 



the Phdf and Phdd polynomial schemes. We can see that the vapor starts increasing 
when the boiling point is reached. The Phdd scheme is as accurate as the Roe scheme 
and captures the correct boiling point, while the Phdf scheme is more diffusive as in 
the saturated boiling channel case and the boiling point obtained by the Phdf scheme is 
inaccurate. 

Without the positivity treatment, the polynomial scheme Phdd also meets some posi- 
tivity problems for vapor volume fractions smaller than = 10^^. The Phdf polynomial 
scheme is more robust due to its diffusivity and allows to compute the case for al = 10^^, 



but the result is not accurate enough to be satisfactory (fig. 9(b)). To obtain an accu- 
rate result even for small void fractions, we therefore use the polynomial scheme Phdd 
additionned with the positivity treatment developed in section |4| called Phdd '■ when a 
positivity problem appears, the step is computed with more numerical diffusion locally 
where the problem is encountered. We can verify on this subcooled boiling channel test- 
case that positivity problems are solved whenever they appear, allowing to compute the 
test-case with very small vapor volume fractions while keeping the result accurate. The 



void fraction profile is displayed on fig. 9(b) for the Ph°^d method with = 10 ^. As the 



diffusion is increased very locally, i.e. only on the faces whose neighbouring cells present 
a lack of positivity, the result remains very accurate and the boiling point is correctly 
captured. The Tank method gives a very good result in terms of stability as it is able to 
compute the test-case with = 10~^ without the positivity treatment and with a very 
good accuracy. However, the computational cost is very high due to the computations of 
the hyperbolic tangent of matrices: the computational time is multiplied by 85 on this 
case compared to the Phdd polynomial scheme. The stability properties of the Tank 
scheme are thus very promising but at the present time it cannot be used in practice due 
to its high computational cost. 

In tab. [2| we provide some statistics on the method Pij^id different initial and 
inlet vapor volume fractions, and for different time-steps of the implicit computation. 
The subcooled boiling channel test-case has been run with CFL=10 and CFL=30, where 
a CFL of 1 corresponds to the stability condition for an explicit scheme. With flmax the 
maximum signal speed, we have: 

At = CFL-^ (5.1) 

I '-'■max I 

First, we can see on tab. [2] that the total number of time-steps where positivity problems 
or difficulties of computation of the pressure have appeared remains small: less than 0.2% 
with CFL=10, and less than 1.8% with CFL=30. The influence of the time-step At on 
the occurrences of positivity problems is very clear, as there is a significant increase of 
problematic time-steps for a CFL of 30 compared to a CFL of 10. At each problematic 
time-step, the method Ph°dd iterate until the right diffusion is found so that the 

scheme is positive. We observe that the average number of iterations is close to 1. This 
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(a) Void fraction - = 10 



(b) Void fraction - q;'„ = 10 ^ 



Figure 9: Boiling channel in the subcooled case with 150 cells. Void fraction for a„ = 10~^ 
(a) and for = 10"^ (b) at time t = 5 s as a function of space. Results by the Phdf 
(green squares) and Prdd (blue diamonds) polynomial schemes, and the Tank (violet 
triangles) method. The Roe scheme (black circles) is depicted on fig. (a) but breaks 
down in case (b). The vertical dashed line indicates the boiling point ybou = 1.21m. 

means that in most of the cases, the positivity problem is solved at the first iteration, i.e. 
with D = 10. More rarely, it takes more than one iteration to obtain the positivity. If 
increasing the diffusion does not solve the positivity problem on a cell, the time-step At 
has to be reduced (by dividing it by 10). We can see that the time-step seldom has to be 
reduced. 



5.3 Tee junction 

The two-dimensional tee-junction test-case shows a dynamic separation between the liquid 
and the vapor phase, thus creating accumulation of vapor in some spots and disappearance 
in others. It consists in a two-dimensional horizontal pipe Ti of length 0.877 m and 
diameter 0.055 m, connected to an other horizontal pipe T2 of diameter 0.055 m in a; = 
0.197 m and whose length from the junction is 0.7196 m. A mixture of water and steam 
enters the first pipe Ti at x = (0.0, 0.0). Due to the difference of density and thus inertia 
between vapor and liquid, most of the hquid continues in the first pipe after the junction 
while a big part of the vapor is deported in the second pipe at the junction. Vapor thus 
accumulates at the junction. The phenomenon is only dynamic as no phase change occurs 
in this test-case. This test-cases allows to test the ability of the scheme to deal with a 
large volume fraction range. 



The model is the two-fluid two-phase flow model presented in section 2.1 The source 



terms included in the case are detailed in appendix A. 3. 2 At the inlet, the following 
values are fixed: m„ = (1.0,0.0) m/s, ue = (1.0,0.0) m/s, = 2650 kJ/kg, he = 1607 
kJ/kg and = 0.45. The pressure is fixed at the outlet. At the outlet of the horizontal 
pipe, p^*'"* = 150 10^ Pa, and at the outlet of the vertical pipe, p^""^* = 149.998 10^ Pa. 
A wall slip boundary condition is prescribed on the walls. The computation is implicit. 
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Vapor volume fraction 
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Table 2: Statistics on the subcooled boiling channel test-case over 5s of computation 
(~ 50000 iterations) for the Phdd scheme 



A coarse mesh with 1149 cells (fig. 10(a)[ ) and a refined mesh with 11315 cells (fig. 11(a) ) 
have been used. 

This case cannot be run with the Roe scheme as positivity problems are met with the 
coarse and the refined meshes. A first improvement is brought by the use of polynomial 
schemes as the Phdd scheme is able to compute the case on the coarse mesh without any 



problem. The result for the void fraction is shown on fig. 10 The result obtained by 
the Phdf scheme is too diffusive to be of interest. However, the Phdd scheme alone is 
not able to compute the case on the refined mesh, due to positivity problems. Therefore 
we use the positivity treatment developed in section |4| All the positivity problems are 
overcome by this method and we are able to show the result of the computation with 



^HDD ^S- l^'^ vapor volume fraction. 

In tab. |3| we provide some statistics on the P^^j^, method for the tee-junction test-case 



on a refined mesh. During the computation, the CFL (eq. (5.1)) increases linearly in a 
lap time of 10 s from CFL=50 to CFL=690. Tab. [3] shows that the number of time-steps 
where positivity problems appear remains very small. The average number of iterations of 
the algorithm remains below two iterations. It means that most of the time the positivity 
problems are solved in one iteration of the P^^^d method, i.e. with a diffusion D = 10. 
Only two time-steps have required a diminution of the time step At in order to obtain 
the positivity. 



6 Conclusion 

In this paper, we have considered numerical schemes for multi-phase flow models when 
one of the phases appears or disappears. The causes of the difficulties that standard 
methods face in this situation have been identified. They are: (i) the loss of hyperbolicity 
of the model when a phase appears or disappears ; (ii) the lack of positivity of the scheme. 
Polynomial schemes have been developed to avoid the use of the eigenvector decomposition 
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(a) Coarse mesh (b) P^^j^, - Void fraction 



Figure 10: Tee junction computed by the Phdd scheme on the coarse mesh. Left: mesh 
used for the computation. Right: vapor volume fraction as a function of space (color 
coded). 




(a) Refined mesh (b) Ph°dd ~ Void fraction 



Figure 11: Tee junction computed by the Phdd scheme on the refined mesh. Left: mesh 
used for the computation. Right: vapor volume fraction as a function of space (color 
coded). 
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Number of problematic time-steps 

(in % of the total number of time-steps) 


0.036 


Average number of iterations 
per problematic time-step 


1.67 


Number of time-steps where the time-step At 
had to be divided by 10 to obtain positivity 


2 



Table 3: Statistics on the tee-junction test-case with a refined mesh and the Ph'dd scheme, 
between Os and 5s (85900 iterations). 

of the Roe matrix and tackle problem (i). A specific positivity treatment has been applied 
to the polynomial solver to treat problem (ii). The resulting method is very robust: large 
ranges of void fraction can now be computed with high accuracy. The method has proved 
effective and accurate on test problems on which standard methods fail. An alternate 
method, based on the hyperbolic tangent function, has also been proposed. It is as 
accurate as the polynomial solver and has shown very good positivity properties without 
requiring any positivity treatment. However, it is computationally too intensive. Future 
work will be concerned with improving the computational cost of the hyperbolic tangent 
method, and combining the polynomial method with the all-speed methodology proposed 
in [9l [H] . The latter will allow to treat situations where some parts of the flow are in the 
small Mach- number regime. 



A Appendix 

A.l Eigenvalues of the two-fluid model 

We investigate the structure of the eigenvalues of the two-fluid system (including the 
energy equations). We recall the method employed in ES]- ^^e finite volume 
framework, the system can be written in the quasilinear form: 

BY ^ ,^,,dV 
— +A„ V — =0. 
ot on 

where n is the normal vector on the considered face. We look for the roots of a polynomial 
P(A), the characteristic polynomial of the An matrix, of degree 2{2 + d), d being the space 
dimension. A straightforward computation leads to the following polynomial 

where Ukn is the the projection on the normal vector of the velocity of phase fc, and P4 
is a polynomial of degree 4. It follows immediately that and uin are some of the 
eigenvalues of the system of multiplicity d. 

For the other eigenvalues, we look for an approximation of the roots of P4 and use a 
perturbation method by introducing the small ratio 

e = — , (A.l) 
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where Um is the projection of the relative velocity on the normal vector and is the 
'characteristic' speed of sound, in the two-phase mixture, given by 



PrnjavPe + aePv] 
Pvpt 



with Cm the mixture sound velocity : c: 



1/2 



2 _ 2 

Pm ' 



2^2 



c:,c: 



7 



avPicj + aip^c^g ' 



aepv + avpe 
aePvCj"^ + 0!vpec~' 



dpk 
dp 



1/2 



and Sk is the entropy. The first order approximation of the two-fluid system eigenvalues 
is 



avpeUyn 


+ aepvUin 


aepv 


+ a^pe 


a^piUvn 


+ aipyUin 


aepv 


+ a^pe 




+ a^PiUin 


aepv 


+ a^pe 


aiPvUvn 


+ a^peuin 


aipv 


+ a^pt 



+ am + 



a^p^ + a^pe 



u^j^a^Pyaipe , 
aepv + a^pe ' 



+ 



1 



aip^ + a^pe 



u';.^a^p^aipi. 
a^Pv + a^pe ' 



(A.2) 



with AP* the interfacial pressure default. The approximate formula of the eigenvalues 
associated with the void waves leads to the hyperbohcity condition 



AP* > 



aip^ + a^pi 

which corresponds to Bestion's model for the interfacial pressure term [4j. Expressions 



(A.2) can be written in the following form 



u„ 



u. 



nap 



-Mr 



-VLr) ■ n 



n-a„ + 0(r). 



Kae a/(5 — l)avaiK 



+ aeK 



a^ + apK 



u,) ■n + 0(e' 



with K = denoting in general a small number. 



A.2 Void fraction and pressure wave eigenvectors of the two- 
fluid model, and asymptotic behaviour 



A first-order approximation in ^ (given by A.l ) of the eigenvectors of the two- fluid model 



has been given in |3l] for a perfect gas of constant 7. Let us recall the expression of the 
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right eigenvectors i?3 and associated to the eigenvalues A3 and A4 in (A.l), and which 
are suspected to collapse when the void fraction tends to zero: 



R: 



3,4 



Pv 
-^3,4 

^3,4 



1 

-I 

7 



1 2^ 

2«J 



Mr 



'2 
P 



:U,, 



A3,4) 



Pv Pi 



(A.3) 



Let us now suppose that the vapor phase disappears and the vapor volume fraction 
tends to zero. In this case, we assume that the relative velocity Um also tends to 
zero. The fast eigenvalues Ai and A2 are now equal to m„ ± am- They remain distinct and 
the eigenvectors associated to these eigenvalues do not collapse. As for the intermediate 
eigenvalues, the void eigenvalues A3 and A4, the form of which are recalled below: 



A 



3,4 



(u£ H u,. ± Mr) ■ n + 0(4 ) 



tend to Un- 

One can also check that the eigenvectors i?3 and i?4 have the following form R 
R^ + 6R + 0{^^), namely: 



R 



7 



1 

_Pl 

Pv 

Un 

pi 

Ur 

Pv 

-u 



-U 



Pv Pi 



± 







^(5 



Urn"' 

^rny P 

Pv 

U ■ Ury/a^(3 





+o{e 



(A.4) 



with /3 



and B — > -i/— when — )■ 0. When tends to zero, tends to 
zero and so does ^. Therefore, bR ~ oihur also tends to zero and -R3 and -R4 collapse. 



A.3 Test-cases 
A.3.1 Boiling channel 

The model used in the boiling channel test-case is the two fluid two phase flow model 



presented in section 2.1 Here are the modeling terms included in the case. We assume 
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that while hi < the hquid saturation enthalpy, the heat flux is only implied in the 
heating of the liquid (heat transfer). When hi > hf^*, the heat flux becomes implied in the 
evaporation only and therefore results in mass transfer. The mass transfer also implies a 
transfer of momentum and energy. All numerical values are indicated below. 



1. The interfacial pressure term is the Bestion's modeling term (2.10) with 5q = 1.1 
and K = 10~^. 

2. Interfacial velocities and enthalpies: 

K = hT\ h\ = hf. 

3. Wall heat transfer concentrations: 



Q"^ = q if h,< hT\ 
= otherwise. 

q: = o. 



4. Mass transfer: 



5. Drag force: 



6. Wall friction: 



7. Gravity: 



r = if hf< h'"' 



i 
L 



i 5 

otherwise. 



F^^ — —F^^ — — -Cz)ajPm|u,.|u,.. 
o 



_ f afcPA:|Ufc|Ufc 

k 



D 



h 



^ext S • 

Numerical data and auxiliary relations are given in tables |4] and [5j 
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Dh = 0.628 m 


Lh = 3.65 m 


NpcH = 10 


uo = 0.7802 m/s 


tti = — - with Tj = 5.10"^ 

ri 


Cd = 0.44 


/ = 0.017 


g = -9.81 m/s^ 





Table 4: Numerical data for the boiling channel test-case 



L = /if* - /if* 


1 1 

Pv Pe 


Ur = Uy — Ui 


Pm = a^pv + aipe 


NpchUqL 

LhViv 





Table 5: Auxihary relations for the boihng channel test-case 



A. 3. 2 Tee Junction 

The model used in the tee-junction test-case is the two-fluid two-phase flow model pre- 
sented in section 12.11 Here are the source terms included in the case: 



1. The interfacial pressure term is the Bestion's modeling term (2.10) with 5q = 1.1 
and K = 10~^. 

2. Interfacial velocity: 

u* = ue, 

3. Drag force: 

o 

with tti = ^,ri = 0.3165 10-^ and Cd = 0.44. 

4. Wall friction: 

_ f ttfcPfc|Ufe|Ufc 



TP" 



with Dh = 1 m and / = 0.05. 
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A. 4 Coefficients of the polynomial Phdf 
The Phdf polynomial is written: 

17 



QkX 



2k 



k=0 



with the Ofe given by: 



ao = 6.209633161688544e - 02 

01 = 4.516480010541272e + 00 
Qi = -3.049057345414379e + 01 

02 = 1.657256844603353e + 02 
04 = -6.133533687894306e + 02 
as = 1.580698142537855e + 03 
og = -2.879210705862515e + 03 
a7 = 3.673105197391366e + 03 
og = -3.121407591514732e + 03 



On 
ai2 
ai3 
ai4 

ai5 
ai6 
an 



1.512887040780976e + 03 
-2.111058506112595e + 02 
9.753698909265717e + 01 
-6.475861637079317e + 02 
8.947647548149256e + 02 
-6.303841204016171e + 02 
2.586951712420909e + 02 
-5.941358894806618e + 01 
5.960406627331660e + 00 
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